A Dipole Polarizable Potential for Reduced and Doped Ce02 from First-Principles. 



Mario Burbano ^, Dario Marrocchelli^Q Bilge Yildiz ^, Harry L TuUer ^, Stefan 
T Norberg Stephen Hull ^, Paul A Madden^, and Graeme W. Watson^ 
^ School of Chemistry and CRANN, Trinity College Dublin, Dublin 2, Ireland 
^ Department of Nuclear Science and Engineering, Massachusetts Institute of Technology 
^ Department of Materials Science and Engineering, Massachusetts Institute of Technology 
Department of Chemical and Biological Engineering, Chalmers University of Technology 
^ The ISIS Facility, Rutherford Appleton Laboratory and 
^ Department of Materials, University of Oxford, 
Parks Road, Oxford 0X1 3PH, United Kingdom 

In this paper we present the parameterization of a new interionic potential for stoichiometric, 
reduced and doped Ce02. We use a dipole-polarizable potential (DIPPIM) and optimize its pa- 
rameters by fitting them to a series of DFT calculations. The resulting potential was tested by 
calculating a series of fundamental properties for Ce02 and by comparing them to experimental val- 
ues. The agreement for all the calculated properties (thermal and chemical expansion coefficients, 
lattice parameters, oxygen migration energies, local crystalline structure and elastic constants) is 
within 10-15% of the experimental one, an accuracy comparable to that of ab initio calculations. 
This result suggests the use of this new potential for reliably predicting atomic-scale properties of 
Ce02 in problems where ab initio calculations are not feasible due to their size-limitations. 



I. INTRODUCTION 

Cerium dioxide, Ce02 or ceria, is an important 
material which has found applications in several tech- 
nologically relevant areas such as catalysii^ and Solid 
Oxide Fuel Cells (SOFCsj^. In catalysis, it plays an 
important role thanks to its oxygen storage capability, 
due to the ready oxidation state change from Ce^+ to 
Ce^^ upon reduction and the reverse upon oxidation^. 
These properties are made use of in Three- Way Catalysts 
(TWC), where the stored oxygen aids in the oxidation 
of CO to CO2 under reducing conditions while, under 
fuel-lean conditions, the reduction of NO to N2 is 
assisted by the uptake of oxygen by ceria. Doping ceria 
with aliovalent cations, such as Gd, Y or La, leads to 
high ionic conductivity in the intermediate temperature 
range (500 - 800° C), thus raising the prospects of 
ceria-based electrolytes for application in SOFCs"* ^ . 

Over the past 5 years, significant progress has been 
made in the description of this material by means 
of ah initio computer simulation^SHUl^ using Density 
Functional Theory (DFT). In particular, the use of the 
DFT-I-U approach, where the U parameter provides an 
improved description of the strongly correlated cerium 
4/ states in partially reduced ceria, has led to a much 
improved understanding of the electronic and structural 
properties of this material. Unfortunately, DFT calcu- 
lations are still severely limited by system size and the 
time-scales that can be studied; this high computational 
cost usually limits this approach to static calculations 
only. For this reason, reliable interatomic potentials 
which allow the study of thousands of atoms on the 
nanosecond scale are desirable. This is particularly 
true for the study of the ionic conductivity of ceria. 
Indeed, the role of grain boundaries in the formation 
of space charge regions, or the vacancy and/or cation 



ordering tendencies, which are responsible for the drop 
in conductivity after a critical vacancy concentration 
(around 3-4 %), are long-range in nature and necessitate 
large simulation boxes. 

In a recent paper, Xu et al^^ compared six differ- 
ent interatomic potentials for ceria available in the 
literatur^i^JtiSI and tested their accuracy by reproducing 
a series of experimental data (lattice constants, thermal 
expansion, chemical expansion, dielectric properties, 
oxygen migration energy and mechanical properties). 
Two main limitations were found. The first was that 
none of the reviewed potentials could reproduce all the 
fundamental properties under study, although some 
displayed higher accuracy than others. While all the 
potentials could reproduce the static properties, such as 
lattice parameters and elastic constants, they all failed 
at reproducing the thermal expansion coefficient, and, 
to a lesser degree, the oxygen migration energy, for pure 
Ce02. Indeed, some potentials gave thermal expansion 
coefficients which were one order of magnitude smaller 
than the experimental one and also severely underes- 
timated the oxygen migration energy. Thermal and 
chemical expansion properties of ceria are particularly 
important in the context of SOFCs given that differential 
expansion of the components has a detri ment al effect on 
the long term durability of the fuel cell^^SEI], A second 
problem evinced from the study by Xu et al. was that 
not all the interatomic potentials have a complete set 
of parameters available for the study of both doped 
and reduced Ce02. The potential by Inaba et aP^, for 
instance, properly reproduces the thermal expansion 
coefficient and the elastic properties, but cannot be 
tested for chemical expansion because it does not have 
parameters for Ce"^"*". 

The first limitation can be easily understood by 
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looking at figure [T] where we show the shape of a typical 
interatomic potential. Such a potential is harmonic in 
the vicinity of the equilibrium position (red curve is a 
parabola fitted to the potential) but at distances away 
from the equilibrium position, it deviates from that 
shape and becomes anharmonic. It is this anharmonicity 
which is responsible for the thermal expansion observed 
in solid materials. A potential's failure to reproduce the 
experimental thermal expansion means that it is not 
properly parametrized and/or the potential's shape is 
not correct at distances greater than the equilibrium 
position. In particular the strong underestimation 
of the thermal expansion coefficient, as observed by 
Xu et al, indicates that the potential maintains a 
harmonic description of the system in regions where this 
approximation is not valid. The potential also deviates 
from the harmonic behavior at distances shorter than 
the equilibrium one, with the interatomic potential 
being more repulsive than its harmonic approximation. 
This feature plays an important role in the calculation 
of the oxygen migration barrier. When an oxygen ion 
hops from one site to another, it has to squeeze between 
neighboring cations, so that the average interatomic 
distance between the oxygen and these cations is much 
smaller than at equilibrium. If a potential's description 
of this interaction is predominantly harmonic, then the 
repulsion between these ions will be underestimated 
and, consequently, the migration energy barrier too will 
be underestimated, as observed by Xu et al^ 
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Figure 1: The shape of a typical interatomic potential (black 
line) and of its harmonic approximation (red line). 

It would seem from the above analysis that the limi- 
tations of the existing potentials are due to the fact that 
they are too harmonic. This might be due to the way 
they are parametrized, given that they are optimized by 
fitting their parameters to a small data set of experi- 
mental properties, usually elastic constants and lattice 
parameters. These are equilibrium or near-to-equlibrium 
properties, which allow the sampling of only a small re- 



gion of the potential energy surface (PES) of the system, 
around the equilibrium distance. Little or no information 
is provided about the long- or short-distance behavior 
where the potential becomes anharmonic. In this paper 
we use a different methodolog;^ which was used success- 
fully for similar oxide systema^SHaD^ iq parametrize inter- 
atomic potentials for stoichiometric, reduced and doped 
Ce02 . The key idea of this methodology is that the po- 
tential's parameters are fitted to a series of DFT calcu- 
lations on high-temperature, distorted Ce02 configura- 
tions. This allows sampling wider regions of the PES 
than those accessed by macroscopic equilibrium observ- 
ables. We obtain parameters for stoichiometric, reduced 
and doped (La, Gd and Y) Ce02 and test them against 
the experimental data. The agreement is quite good for 
all the studied properties, including thermal expansion 
and oxygen migration energy barriers. This methodol- 
ogy can be easily extended to other dopant cations in 
ceria or similar materials. 



II. POTENTIAL DEVELOPMENT 

The interaction model used in this work is the same 
as that used in previous work on si milar systems, such 
as GeOjsnmi^ ^oped ZrOspTMUllIl^ and LiaCP. The 
model ^, known as DIPole Polarizable Ion Model (DIP- 
PIM), includes a pair potential (a Buckingham term plus 
Coulombic interactions), together with an account of 
the polarization effects that result from the induction of 
dipoles on the ions. This model is conceptually similar 
to the shell model used by many authorsSHUl Here we 
use formal ionic charges (0^~, Ce^"*", Ce^+, La'^+, Gd"^"*") 
which should ensure better transferability. A description 
of this model and the notation used for its parameters 
is reported in Appendix |^ The parameters for these 
potentials were obtained by matching the forces/dipoles 
obtained from DIPPIM to first-principles reference 
datcP3. Such an approach has been applied successfully 
in the case of other oxide material^^UM] jj^ ^^le following 
we give a brief description of the first-principles-based 
reference calculations and the force/dipole-matching 
procedure. 

First, a potential was parametrized for Y doped 
Ce02, since this system will be the object of a detailed 
experimental and computational study by the authors 
of this papei"^^ . This was done by performing DFT 
calculations, in the local density approximation (LDA), 
using the CPMD cod^^. For this system, the simple 
DFT approach is known to give the correct valence states 
(Ce''+ and Y^+), so that a DFT-hU functional was not 
needed. Twelve 2x2x2 supercells with Ceo.5Yo.5O1. 75 



This model is implemented in an in-house molecular dynamics 
code called PIMAIM. 
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compositions and a total of 88 atoms were constructed. 
Each model supercell was obtained from high tem- 
perature (2500 K) MD simulations that were run for 
50 ps in order to reach structural equilibrium. The 
forces on each species were determined directly from 
each DFT calculation, and the dipoles were obtained 
from a Wannier analysis of the Kohn-Sham (KS) wave 
functions'^^. Once the information about both forces and 
dipoles was gathered from the ab initio calculations, 
the parameters in the interatomic potential were fitted 
to them. One problem with DFT calculations is the 
uncontrolled representation of the dispersion interaction. 
Although dispersion energies constitute only a small 
fraction of the total energy, they have a considerable in- 
fluence on transition pressures and, in particular, on the 
material density and stress tensor. For this reason, the 
dispersion parameters were not included in the fits but 
were added afterwards, as discussed by Madden et al^. 
The parameters from refs; ^-' ^ ^' were used. The resulting 
parameters for the potential are reported in table [ij The 
short-range parameters for the 0-Y interaction are in 
line with those obtained in a previous studjESl and the 
value for the 0^~ polarizability is close to that obtained 
from independent ab-initio calculation^^. 

A potential was parametrized for reduced Ce02 as 
well. Similarly as above, first a short ab initio MD 
simulation was performed on a 2 x 2 x 2 supercell 
with CeOi.875 composition at high temperature. This 
configuration was used to calculate the forces acting on 
the ions and these, in turn, were used for the fitting 
procedure. These calculations were performed with 
the Vienna Ab-initio Simulation Package (VASP^PSD 
within the DFT-I-U framework. A value of U = 7 eV was 
chosen to ensure correct localization of the / electrons. 
Although this is slightly higher than is generally usecP^l, 
it was found to be necessary to ensure localization in 
the non-equilibrium structures used in the potential 
fitting. We used the Generalized Gradient Approx- 
imation (GGA) with the Perdew-Wang 91 (PW91) 
exchange-correlation functional and an energy cut-off 
of 400 eV. The Ce'^'*' cations were identified as those 
with spin 1, while the Ce^^ ones have no spin. The 
Ce^+ and Ce*"*" cations are then treated as two different 
cationic species in the potential parametrization and 
separate terms are obtained to describe the interactions 
between, for instance, Ce'*+ - and Ce^+ - O^^. 

The parameters for the Ce'*+ - O^^ and O^^ - O^^ 
were then fixed to the values obtained before for the 
Y doped Ce02 system, so that these parameters are 
consistent and can be used all together. This procedure 
is equivalent to treating the Ce^+ cations as a dopant 
species, in the same way as was done for Y. While an 
ionic model for reduced Ce02 might not be justified a 
priori^ the picture arising from DFT-I-U calculations 
and experiments seems to confirm this model. The / 
electrons are, indeed, found to be strongly localized on 
the Ce'^"'" cations in agreement with an ionic picture. 



Also, a recent study on Ce'^"'"/Ce^+ ordering in ceria 
nanoparticle^^] showed that strikingly similar relative 
energy ordering of the isomers and atomic scale struc- 
tural trends (e.g., cation-cation distances) are obtained 
in both the DFT and interionic-potential calculations, 
which, again, proves the validity of this approach. The 
resulting parameters for the Ce^"*" cations are reported 
in table [l] The potential parametrization was performed 
using two DFT codes because the CPMD package does 
not have an implementation of the DFT-I-U framework 
needed to describe reduced Ce02, while the functionality 
required for the Wannier analysis was not available in 
VASP at the time. 

Finally, the same procedure can be repeated to include 
other cations. In this case we obtained parameters for 
La'^"'" and Gd'^^. This was again done by performing a 
short MD simulation on Gd and La-doped Ce02 at high 
temperature and using the final configuration to calculate 
forces to which the potential parameters can be fitted. 
Once again, we fixed the parameters for the Ce*+ - 0^~ 
and 0^~ - 0^~ to the values obtained before for the Y 
doped Ge02 system. This ensures that these parameters 
are all consistent and allows the study of mixed systems 
- such as, for instance, partially reduced, and Y and La 
doped Ce02. The resulting parameters for the dopant 
cations are reported in table [T] This procedure can be 
repeated for as many dopant cations as desired. 



III. RESULTS 

The quality of the potentials can be assessed by com- 
paring model predictions to experimental data, since no 
experimental data was used in the optimization of the 
model parameters. In this section we present our predic- 
tions, using the herein developed interatomic potentials, 
on lattice parameters, local crystalline structure, thermal 
and chemical expansion, oxygen migration energies and 
elastic constants, and their comparison with the exper- 
imental values. The way these quantities are calculated 
is described in each section. 



A. Thermal expansion 

Xu et al}^ find thermal expansion to be one of the 
most difficult properties to model in pure ceria. Figure [2] 
shows the calculated expansion of the lattice parameter, 
in the 300 - 1000 K temperature range. This is defined 
as (a-ao)/ao, where a is the lattice parameter for a cer- 
tain temperature and ao is the lattice parameter at 300 
K. This was obtained by performing molecular dynamics 
simulations with 4x4x4 supercells in an NPT ensem- 
ble, at the required temperatures. We used baros tats 
and thermostats as described by Martyna et al^^^ and 
we set the external pressure to zero. The lattice param- 
eters were averaged over a 0.1 ns long simulation and 
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Table I: Parameters of the DIPPIM potential. All values are in atomic units. A description of the parameters is reported in 
[A] For the polarizability part of the potential we report only those parameters with b not equal to zero. 





O-O 


Y-0 


Ce-'+.O 


Ce^+-0 


Gd-0 


La-O 




55.3 
6.78 
50000 
0.85 


111.1 
1.377 
50000 
1.35 


105.9 
1.269 
50000 
1.4 


218.7 
1.473 
50000 
1.35 


236.9 
1.566 
50000 
1.35 


98.6 
1.257 
50000 
1.35 


CI' 

disp 


53 
1023 
1.0 


12 
240 
1.5 


12 
240 
1.5 


12 
240 
1.5 


12 
240 
1.5 


12 
240 
1.4 



ao2- 14.9 

ay3+ 2.60 

ace*+ 5-0 

ace3+ 11-2 

«Gd3+ 6.8 

ai„3+ 10.8 



^o-o 
o-o 

,0-Ce3 + 

"d 

y3+_y3 + 
-i^3+_v3+ 



'd 

,La^+-0 
■D 



1.73 
0.45 
1.82 
2.92 
1.67 
0.90 
1.93 
2.81 
1.69 
-1.21 



0-Ce'^ + 
D 

0-Ce* + 
D 

Ce^+-0 
D 

Ce^+-0 
D 

D 



y3+. 



o 
-o 
-o 



1.76 
1.93 
1.82 
-2.50 
0.59 
-.33 
1.94 
-.44 



Ce*+-0 
D 

Ce^+-0 
D 

0-y3+ 
D 

0-y3+ 
D 

D 

y3 + 
D 

6°- 



3 + 



— La' 

D 

iO-La3 + 



1.76 
-.47 
1.67 
1.62 
0.75 
-.45 
1.69 
2.02 



used to calculate the percentage expansion. The corre- 
sponding thermal expansion coefficient is extracted by 
fitting a straight line to our data (see figure [2|. The ob- 
tained value is a = 9.0 x 10^® K^^ which is within 20 % 
of the experimental values (see table |ll]). This is a sub- 
stantial improvement compared to the GrimeJ^^, Gotte 
2004^^ and Gotte 2007*131 potentials which gave a ther- 
mal expansion coefficient of 1.27 x 10~^ , 6.65 x 10~^ 
, 7.31 X 10^^ K^^, respectively. We remind the reader 
that no empirical data was used in the parameterization 
of this potential, so that this potential was not manually 
optimized in order to reproduce the thermal expansion 
coefficient. 



Thermal expansion coefficient (10 ^ K ^) Reference 



9.0 


This work 


10.7 


Ref.133 


11.1 


Ref.Sa 


11.6 


Ref.E3 



Table II: Comparison between the experimental and simu- 
lated thermal expansion coefficients. 



B. Elastic properties 

Elastic constants and the bulk modulus were calcu- 
lated and compared to experimental values. The three 
independent elastic constants were obtained by straining 
an optimized simulation cell in different directions by 
a small amount (typically a fraction of a percent) and 
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Temperature (K) 



Figure 2: Simulated lattice expansion (points) versus temper- 
ature for stoichiometric Ce02 . The red line represents a linear 
fit of the MD data. The thermal expansion, a, extrapolated 
from the fit is 9.0 x 10"® K'V 



which can be obtained from the elastic constant^^. Fig- 
ure |3] reports the Young's modulus for reduced ceria as 
a function of the lattice constant, together wit h the ex- 
perimental data from the work of Wang et al^^^ and 
the ones obtained with the Grimes and Gotte 2007 po- 
tentials, as reported by Xu et al^. The scale of the plot 
is the same as in Fig. 3 in refl^ to facilitate a compari- 
son. Our interatomic potential shows the expected elastic 
softening as a function of the lattice constant (or analo- 
gously, as a function of non-stoichiometry) . When com- 
pared to the experimental data, it seems that our poten- 
tial slightly overestimates this softening, although most 
of the values obtained with our simulations are within 
experimental error. The level of accuracy is indeed bet- 
ter than the Grimes potential and comparable with the 
Gotte 2007 potential. Finally, preliminary calculations 
on Cei_^M^02-:r/2, with M = , Gd^+ and La3+, 
indicate that this system shows a similar softening as a 
function of the dopant concentration, x. 



the resulting stress tensor is recorded after relaxation 
of the atomic positions. After repeating this procedure 
for several magnitudes of positive or negative strain, the 
linear relationship between strain and stress is used to 
obtain the elastic constants at K. The bulk modulus 
was extracted from a volume versus pressure curve. 



Table |III| reports the three elastic constants and the 
bulk modulus. These are within 15% the reported exper- 
imental values, with the exception of Cn, whose value is 
overestimated by 32%. We believe this agreement to be 
quite good, especially considering that the experimental 
data for these properties are quite scattered. 




Lattice Constant (A) 



Property (GPa) MD ExperimentaPEll 



Cii 


552 


386 


- 450 


Cl2 


137 


105 


- 124 


C44 


66 


60 


- 73 


B 


275 


204 


- 236 



Table III: Experimental and simulated elastic constants and 
bulk modulus for stoichiometric Ce02. 

One of the main motivations of the work of Xu et al^ 
was to find a reliable potential to describe - and per- 
haps explain - the observed elastic softening of ceria with 
decreasing oxygen partial pressure {Po2)- We therefore 
calculated the Young's modulus by using the following 
formulcPSl, 

3B + G' 

where B is the bulk modulus and G the shear modulus. 



Figure 3: Young's modulus versus the lattice parameter in 
Ce02-x- Red symbols are the experimental data from Wang 
et alW^, green and blue symbols are the simulated data ob- 
tained with the Grime^^ and Gotte 200'i^ potentials, while 
the black ones are from our own potential. 



C. Structural properties of reduced Ce02 

The local crystal structure of reduced ceria has been 
recently studied by neutron diffraction''^^. From that set 
of data, total radial distribution functions were extracted 
for different values of x in Ce02-2;. Total radial distribu- 
tion functions can be expressed in terms of the individual 
partial radial distribution functions, gij(r), weighted by 
the concentrations of the two species, Ci and Cj, and their 
coherent bound neutron scattering lengths, hi and 6j, so 
that 
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n n 

where n is the number of ionic species. The partial 
radial distribution functions gij{r) are given by: 

with nij{r) equal to the number of atoms of type j 
located at a distance between r and r + Ar from an atom 
of type i and pj is the number density of atoms of type 
j, given by pj = CjPQ. These partial radial distribution 
functions can be easily calculated from the simulation 
output. 

In figure |4j we report the calculated total radial distri- 
bution functions, G(r), for different non-stoichiometries 
and we compare them with those extracted from the neu- 
tron diffraction data. The agreement is good for all the 
studied compositions. A visual analysis of the oxygen- 
oxygen radial distribution functions (not shown) shows 
an increased broadening as a function of x, which is in- 
dicative of an increased disorder within the anion sublat- 
tice, in agreement with experiments^. 



O 3 





A 




V 


A- 




CeO 
























4 




6 


8 


10 


12 








V 
















5 






10 






1 










1 






2 


4 




6 


8 


10 


12 




1 










1 






2 


4 




6 


8 


10 


12 



r(A) 



Figure 4: Radial distribution functions, G{r), for different 
values of x in Ce02-a;, at 1273 K. Red lines and black empty 
dots correspond to the MD and experimentaP^l q^^^s respec- 
tively. 



strain in the cell and can eventually cause fracture. For 
this reason, in this section, we test the ability of our 
potential to describe this behavior accurately. In figure [3] 
we report the calculated lattice parameter as a function 
of the oxygen non-stoichiometry in Ce02 at 1273 K and 
compare this with the neutron data from Hull et alJ'^ at 
the same temperature. The agreement is excellent and 
the simulated chemical expansion coefficient, 0.338 A, is 
within 7 % of the experimentally determined chemical 
expansion coefficient, 0.362 A. Such a good agreement 
is encouraging and also indirectly confirms that our 
"ionic" approach, in which we see Ce^+ as a different 
cation species, carries the correct physics. 



▲ MD 




0.05 0.1 0.15 0.2 0.25 

X in CeO, 

2-x 



Figure 5: MD lattice parameter (red triangles) as a function 
of the nonstoichiometry, x, in Ce02-x, with the correspond- 
ing linear fit (red line). The simulations were performed at 
1273 K, the same temperature as the reported experimental 
datJ^ (black circles). The black line is a linear fit to the 
experimental data. 



In table |IV| we report the calculated lattice parame- 
ters for Ceo. 8-^^0.20 1.9 at room temperature, where M 
= Gd'^"'", La'^^ and compare them with the experimental 
values. The agreement is within 1% and the trend of in- 
creasing lattice parameter with increasing cation radius 
is properly reproduced. 



Lattice parameter MD (A) Experimental (Ajp' 



D. Chemical expansion 

The conditions found at the anode side of SOFCs lead 
to the reduction of Ce^+ to Ce^"*" with a subsequent 
change in the lattice parameter. This chemical expansion 
affects the performance of the electrolyte as it creates a 



Gd'^+ 5.426 5.423 

La^+ 5.494 5.476 



Table IV: Comparison between experimentaP^ and simulated 
lattice parameters for Ceo gAfo 2O1 9 at room temperature (M 
= Gd='+, La^+). 
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E. Oxygen migration energies 



In pure ceria, the activation energy (Ea) which 
determines the ionic conductivity at a given temperature 
is composed of the vacancy formation enerergy (Ef) and 
the oxygen migration energy {Em)- In this case, the low 
concentration of oxygen vacancies results in a low ionic 
conductivity. Aliovalent doping generates vacancies 
through charge compensation so that the activation 
energy then becomes the sum of the the migration en- 
ergy plus the association (binding) energy {E},) between 
the dopant cations and the vacancies .'531E2l^ Molecular 
dynamics simulations on the Cei_2:Y2.02_a;/2 system 
were performed for different values of x in order to 
evaluate the diffusion coefficient and the relative oxygen 
migration energies. These simulations were performed 
on 4 X 4 X 4 supercells at a various Y2O3 concentrations. 
Each concentration was initially equilibrated at a tem- 
perature of 2073 K for 120 ps, using a timestep of 1 fs. 
The diffusion coefficients were then extracted from the 
oxide ion mean squared displacement for temperatures 
of 1073 K and above, as explained by Norberg et al.^^. 
The diffusion coefficients were then plotted as a function 
of the temperature and the oxygen migration activation 
energy was extracted. In principle, this activation 
energy contains both the migration energy and the 
association energy between the vacancy and the dopant 
cation. However, in doped-ceria conductors at high tem- 
peratures, which is the situation of interest here, isolated 
vacancies migrate freely so that the activation energy 
is equal to the migration energy only^^^H^. Figure |6] 
therefore shows the predicted activation energy (oxygen 
migration energy) as a function of the Y concentration, 
X. The observed behavior of an increasing migration 
energy as a function of Y concentration is in excellent 
agreement with the experimental findings of Tian et al.^^. 

If the migration energy in figure [6] is extrapolated to x 
= 0, a value of approximately 0.47 eV is obtained. This 
value corresponds to the oxygen migration energy at in- 
finite dilution, i.e. when the oxygen migration energy is 
unaffected by the dopant cations, and compares well with 
the experimental results and with previous DFT calcula- 
tions, as shown in table |Vj In conclusion, this potential 
can successfully reproduce the ionic conduction proper- 
ties of this material. A more detailed analysis of the ionic 
conduction mechanisms in doped Ce02 and the factors 
that affect it will be the subject of a subsequent papei^^. 



> 



Figure 6: Activation energy (black points) in Cei_a;Ya:02_i:/2 
versus Y concentration, x, with linear fit (black line). These 
energies were extracted from an Arrhenius plot of the MD 
diffusion coefficients at high temperatures. 



Oxygen migration energy (eV) Reference 


0.47 


This work 


0.52 


Ref.^^ 


0.40 




0.47 


Ref.[^ 



Table V: Comparison between experimental and simulated 
oxygen migration energies for Ce02-a;, for a; — >• 0. 



We used a dipole-polarizable (DIPPIM) potential and 
optimized its parameters by fitting them to a series of 
DFT calculations. The resulting potential was tested 
by calculating a series of fundamental properties for 
Ce02 and by comparing them to experimental values. 
The agreement for all these properties (thermal and 
chemical expansion coefficients, lattice parameters, 
oxygen migration energies, local crystalline structure 
and elastic constants) is very good, with the calculated 
values being generally within 10-15% of the experimental 
ones. We note that such accuracy is comparable to 
that of DFT calculations, but the computational cost is 
reduced significantly. 



IV. CONCLUSIONS 

Driven by the recent work of Xu et ai, which showed 
that none of the potentials for Ce02 reported in the 
literature could reproduce all the fundamental proper- 
ties of this system, in this paper we have presented the 
parameterization and the accuracy of a new interionic 
potential for stoichiometric, reduced and doped Ce02. 



These potentials can be used to predict and elucidate 
the atomic-scale properties of Ce02 in situations where 
DFT calculations are not practical due to their size lim- 
itations. With such potentials, nanosecond long simu- 
lations on 1000s of atoms can be performed and these 
can be used to understand the structural, chemical, me- 
chanical and conducting properties of this material, as 
previously done for similar system^^SHUl, 
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Appendix A: The DIPPIM model 

In this section we report a brief description of the po- 
tential model em ployed in this work. The reader is re- 
ferred tcp2l23ISll and references therein for further read- 
ing. The interatomic potential is constructed from four 
components: charge-charge, dispersion, overlap repulsion 
and polarization. The first three components are purely 
pairwise additive: 

ym^Y^m. (Ai) 

where qi is the formal charge on ion i. The dispersion 
interactions include dipole-dipole and dipole-quadrupole 
terms 

Here Cq and are the dipole-dipole and dipole- 
quadrupole dispersion coefficients, respectively, and the 
fl^ are the Tang-Tonnies dispersion damping function, 
which describe short-range corrections to the asymptotic 
dispersion term. The short range repulsive terms are 
approximately exponential in the region of physical inte- 
rionic separations. The full expression used here for the 
short range repulsion is: 

y^^P = ^ — + ^ S'^e-'"'-'^ , (A3) 

2<j ^-^ i<j 

where the second term is a Gaussian which acts as a 
steep repulsive wall and accounts for the anion hard core; 
these extra terms are used in cases where the ions are 



strongly polarized to avoid instability problems at very 
small anion-cation separation^^l^ xhe polarization part 
of the potential incorporates dipolar effects only. This 
reads: 

(A4) 

T(i)(r) = -rjr' T^J^ (r) = (3r„r^ - r'6^f,)/r'. 

(A5) 

The instantaneous values of these moments are obtained 
by minimization of this expression with respect to the 
dipoles of all ions at each MD timestep. This ensures 
that we regain the condition that the dipole induced by 
an electrical field E is aE and that the dipole values are 
mutually consistent. The short-range induction effects on 
the dipoles are taken into account by the Tang-Toennies 
damping functions: 

fl^ir^i) ^ 1 - c^^e-^"-^- ^ i^^. (A6) 

The parameters determine the range at which the 
overlap of the charge densities affects the induced dipoles, 
the parameters determine the strength of the ion re- 
sponse to this effect. 
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